Shape Instabilities in the Dynamics of a Two-component Fluid Membrane 



u , 
Oh! 

< 
(N 



O 



o 
o 



o 

On 



I 

o 
o 



X 



p. B. Sunil Kumar and Madan Rao 
Institute of Mathematical Sciences, Taramani, Chennai (Madras) 600 113, India 

(February 1, 2008) 

We study the shape dynamics of a two-component fluid membrane, using a dynamical triangulation 
monte carlo simulation and a Langevin description. Phase separation induces morphology changes 
depending on the lateral mobility of the lipids. When the mobility is large, the familiar labyrinthine 
spinodal pattern is linearly unstable to undulation fluctuations and breaks up into buds, which move 
towards each other and merge. For low mobilities, the membrane responds elastically at short times, 
preferring to buckle locally, resulting in a crinkled surface. 
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A mixture of phospholipid molecules, e.g. DMPC -I- 
SOPC, aggregate in an aqueous solvent to form closed 
bilayered fluid membranes of dimension 10 — 20 /zm. 
Such mixed vesicles exhibit definite equilibrium shapes 
depending on the ambient temperature, osmotic pres- 
sure, relative concentration and history of preparation. 
Amongst the multitude of equilibrium shape transitions 
predicted by theory and observed in phase contrast 
video micrography, the most dramatic is the phase seg- 
regation induced budding ||l| at which a large 'parent' 
vesicle sprouts out a spherical bud consisting of one of 
the lipid species, attached to it by a narrow umbilical. 
Such equilibrium budding has been oberved in a mixture 
of natural brain (sphingomyelin) lipids jsj . 

Apart from these equilibrium studies, there has been 
relatively little work on the dynamics Q of mixed fluid 
membranes . In the context of a 2-componcnt fluid 
membrane undergoing phase separation, the study of 
nonequilibrium shapes becomes essential, since the lat- 
eral diffusion coefficient of lipids 10~* — lO^^cm^/sec, 
are such that the time scale over which local lipid concen- 
tration relaxes is comparable to shape relaxation times. 
In this Letter we study the shape dynamics of a 2- 
component, open fluid membrane experiencing phase sep- 
aration, using both a dynamical triangulation monte 
carlo simulation (DTMC) and a Langevin description. 

The shape dynamics of a 2-component fluid membrane 
involves 6 slow variables, the conserved total lipid den- 
sity p — PA + Pb, the conserved relative concentra- 
tion (j) = (pa — Pb)/p, the 'broken symmetry' shape 
variable R(iti,M2), and the total momentum density tt, 
which involves both lipid and solvent hydrodynamics. 
The primary dissipation mechanisms arise from the in- 
plane lipid viscosity and out-of-plane solvent viscosity . 
We include the effects of solvent hydrodynamics in our 
Langevin analysis, but not in our DTMC. 

Over length scales larger than the bilayer thickness 
{d ~ 40A), a 2-component fluid membrane can be de- 
scribed by a continuum hamiltonian |]l],^ Ti. ~ Tic + 
Ti.c-tfi + 'H^ + Ti.p, where the curvature energy Q, written 
in terms of the extrinsic curvature H reads {g = det{gij), 
where gij is the metric) 



Tic + 'He 
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while Titj, takes the usual Landau-Ginzburg form, 
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The spontaneous curvature IIq{4>) — co(l + 4>)/'^i re- 
flects the shape asymmetry between the two lipid species, 
and biases the local curvature to be cq or 0. The up- 
per length scale cutoff is set by the persistance length 

« (iexp(47rKc/3fcBT), beyond which thermal fluctua- 
tions drive the membrane into a self-avoiding branched 
polymer phase (kc ~ lOOfcsTroom for DMPC). In prin- 
ciple Kc depends on the local concentration 0, but we 
shall ignore this detail (which allows us to drop the in- 
trinsic curvature term usually present in Eq. The 
energy cost for density fluctuations 5p along the mem- 
brane is Hp = Xo^ J{Sp/2)'^y/g d'^u (where Xo/Po is the 
compressibility of the membrane). 

Our model 2-component membrane consists of a 
fixed number N of two types of hard beads A and B (ver- 
tices), each of diameter a, with N = Na + A^b, linked 
together by flexible tethers so as to triangulate a 2-dim 
open surface embedded in Ti? . The length of each tether 
lies between a = Imin < I < Imax = \/3i (imposes self- 
avoidance locally 1^ ) . The coordination number of every 
bead is restricted to lie between 3 and 9. The configura- 
tions of our model 2-component membrane are weighted 
by the discrete form of Ti. The discretised curvature en- 
ergy Eq. |l| is written as 
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where the index i denotes a vertex, (ij) the tether con- 
necting i and j and sum over all tethers ema- 
nating from i. The term = 1 — • , where ha 
and are the local outward normals to the triangles a 
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Fig. 1 Configuration snapshots from the DTMC simulation at times t = 4000 , 6500 , 17000 in units of MCS, for a critical 
quench. The number of particles is A'^ = 1140, and the number of fiips is Nf = ISA'^. The labyrinthine pattern, typical of 
spinodal decomposition, break up into buds at late times. 



and P sharing a common tether (ij). The intrinsic cur- 
vature Ki, calculated from the deficit angle at i, must be 
added to -ff(y) to obtain Eq. |^ in the continuum limit. 
The ± sign in the third term, denotes the sign of the local 
mean curvature — it is positive if the outward normals to 
adjacent triangles point away from each other and neg- 
ative otherwise. The concentration (pi takes values ±1 
depending on whether the vertex i is occupied by A or 
B. In the second term of Eq. |[ is the sum of the areas 
of the triangles with i as the vertex. The discrete version 
of Eq. H is clearly. 
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where the sum < ij > is over vertex pairs connected by a 
tether. Finally, Tip is contained in our imposition of the 
hard sphere and tether constraints. 

The DTMC for a 2-component fluid membrane consists 
of bead movements in TZ^ to change the shape ||] , tether 
flips to simulate fluidity (dynamical triangulation) 
and exchange moves between A and B to promote phase 
segregation ||] . Moves which violate the hard sphere and 
tether constraints, graph planarity (no tether intersec- 
tions) and topology of the surface |9| are rejected. The 
bead moves and the tether flips are made according to 
a Metropolis algorithm with the hamiltonian Eqs. |3[ ^. 
The exchange moves between A and B are implemented 
by the usual Kawasaki dynamics [lO[ |, with a transition 
probability W{i ^ j) = [1 - tanh {AE/2kBT)]/2, where 
AE is the energy difference between the final and the 
initial configuration. One monte carlo sweep (MCS) is 
defined as iV attempted bead moves, during which we 
make Nf attempts at flipping tethers and Nex attempts 
at Kawasaki exchanges. We shall work with two differ- 
ent boundary conditions — rigid boundary conditions 



on a hexagonal frame (fixes the projected area) and free 
boundary conditions (unframed), in either case the teth- 
ers along the boundary are not flipped. Details of this 
algorithm may be found in 

A variant of this algorithm was used to study the phase 
separation dynamics of dense binary fluids confined to 
two dimensions ||Tl[| . The growth laws and the form of 
the correlation functions were consistent with a variety of 
other simulations jl^ . We note that the lateral mobility 
of the beads may be tuned by Nf. 

At time t = 0, we start with a homogeneous mixture 
(Na ^ Nb^ N/2 ; N = 1140) of A and B hpids on a flat, 
randomly triangulated 2-dim surface. After a quench into 
the coexistence regime [ p^ , the fluid membrane evolves 
by the monte carlo dynamics just discussed. We shall 
investigate the generic case when the time scales, r and 
T0, over which the shape and the concentration change 
are comparable {Nf,x = N), relegating a more detailed 
analysis to later ||l^ . The values of parameters have been 
chosen to be J = Kc = 'to = 1 in units of fcsT = 0.25. 

Before showing the sequence of monte carlo configura- 
tions, let us discuss the very early time spinodal instabil- 
ity, when we are justified in using a Monge parametriza- 
tion of the surface 'R{ui,U2',t) = {x,y,h{x,y;t)). The 
continuum Langevin equations of motion for h and 
read. 
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We perform a linear stability analysis about a flat mem- 
brane with h{x,t) = ^q^^^i^)^^'^^ ^^'^ 4>{^i't) = 00 + 
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X]q S4>q{t)e^'^'^ . The overdamped modes are unstable for 

q < q*, where q* = - 30q + cI/kc)/cf- The small 
q dispersion of the unstable mode is easily worked out 
— Lo ^ q^ for 0g < 1/3 and uj ^ q^ for 0g > 1/3. Sol- 
vent hydrodynamics, introduced through a rcnormalized 
T ~ {vq)~^ iv = solvent kinematic viscosity), changes 
the dispersion to cj ^ q^, leaving the other features 
unaffected. The coupling to shape fluctuations shifts the 
(mean field) spinodal line closer to the coexistence line : 
(j)Q — 1/3 + Cq/3kc. The configuration coarsens into the 
usual labyrinthine pattern. Fig. 1(a), which is a snapshot 
from our DTMC at t = 4000 MCS. 

Subsequent evolution depends on the magnitude of the 
lateral mobility. Let us first consider the case when 
Nf = 15N, which corresponds to mobile lipids in the 
fluid phase. A look at a time sequence of configurations 
generated by our simulation Fig. 1 , reveals an unusual 
scene. The bicontinuous patches, which are a generic fea- 
ture of binary fluids at a critical quench, break up into 
disconnected buds as time progresses ! Further coars- 
ening occurs through the motion of these buds, as they 
approach each other and subsequently coalesce. The late 
time behaviour depends on the final equilibrium config- 
uration of the membrane. 

To understand the breakup of the bicontinuous patches 
(strip instability) , we again work within the Monge gauge 
and perform a linear stability analysis about a flat strip 
of width 2uo with 0=1 within the strip and = — 1 
without. We parametrise the perturbed configuration by 
h ^ ha + ^q(5/iq(t)e*(9i^+922') and (j) = tanh[(M(x,t) - 
y)/^] + tanh[(u(a;, t) + y)/^] — 1, where u{x,t) = uq + 
S<ji '^^91 (Oe*'^"^ and ^ = y/2a is the interfacial thickness. 
The time evolution of the fourier amplitudes evaluated at 
y — ±uo is, to linear order, 
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Such a perturbation leads to a long wavelength insta- 
bility of the strip. The dispersion of the unstable mode 
is uj ql, a,s qi 0. Hydrodynamic coupling to the 
solvent via a rcnormalized r(q) does not alter this disper- 
sion. This instability, arising from a competition between 
the interfacial and the curvature energy |]l| , leads to local 
outward budding since co > 0. This linear instability will 
grow into the complete buds seen in Fig. 1(c), connected 
to the rest of the membrane through a narrow neck. The 
typical size of the bud is Rh — Hc/ {cq + a). The cq = 
case is special since the linear analysis does not predict a 
strip instability ; nonlinear terms however would initiate 
an instability to budding. Since the bud could form ei- 
ther outward or inward, the dynamics of bud formation is 
logarithmically slow as can be seen by a mapping onto a 
1-dim, scalar time-dependent Ginzburg-Landau (TDGL) 
dynamics JTof . 




Fig. 2 Configuration snapshot from the DTMC simula- 
tion at t = 17000 MCS for a critical quench. The number of 
particles is the same as in Fig. 1, but Nf — N. The late time 
configuration is now 'crinkled'. 

When the mobility is small {Nf = N), the interme- 
diate time shapes are dramatically different. Instead of 
erupting buds, the membrane develops sharp ridges at 
the A|B interface and appears crinkled (Fig. 2) . We ar- 
gue that these crinkles are a consequence of a disclination 
induced buckling instability. The A|B interfacial energy 
can be reduced both by shrinking the interface perime- 
ter and by lowering the local coordination number of the 
beads at the interface. This introduces a preponderance 
of positive orientational disclination defects at the inter- 
face, in an otherwise 6-fold coordinated network. Though 
the membrane is still in the fluid phase jist , the high lipid 
viscosity rji results in a nonzero shear modulus fi ^ Vi/^s 
over time scales smaller than the shear relaxation time 
p6[ . The membrane therefore responds elastically to the 
presence of these disclinations. Placing a positive discli- 
nation on a flat elastic sheet of size R costs an energy 
proportional to KqR^ , where Kq = A:Ix{ijl -\- X)/{2fi + A) 
is related to the elastic moduli. Buckling eliminates the 
elastic contributions to the stress field created by this de- 
fect |l^, leaving only a bending energy cost which goes 
as Kc In R. As shown in |p7[] , the elastic membrane buck- 
les when KqR^ > 160 Kc. This physics is consistent with 
our monte carlo shape evolution of a flat membrane hav- 
ing a single patch of A in a sea of B lipids with Nf = N. 
The membrane goes into a buckled conical shape, with 
the A species at the tip. This buckled shape is clearly a 
long-lived nonequilibrium effect. The same elastic con- 
siderations would imply that negative disclinations would 
induce a long-lived negatively curved (saddle) buckled 
surface. The relative fraction of positive and negative 
disclinations depends on the boundary conditions at the 
frame. Unless there is an excess charge, either due to 
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initial conditions or by having an open boundary, the 
membrane will be asymptotically flat. The crinkled look 
of Fig. 2 is a consequence of a random distribution of 
such positively and negatively buckled surfaces. 

The dynamics of coarsening and shape changes can 
be studied from the behaviour of a variety of two- 
point correlators. We define as Re, the length 
scale associated with the interfacial energy density, 
{E) ^ N^'^Y.<i.j>(^ ~ 4'i4>3), and Rh, the length ex- 
tracted from the total mean curvature density (H) = 



In addition, we mea- 



sure R(j,(j, from the first zero of the intrinsic equal- 
time correlator for the local concentration Ti{r2,t) = 
W Su i't'i^ + r2, (/"(u, t)) and Rhh from the first zero 
of the intrinsic equal-time correlator for the local curva- 
ture A,(r2,i) = J2uiH{u + r2,t)H{u,t)), where \r2\ 
is the geodesic distance between points u and u + r2 de- 
fined on the 2-dim manifold. Geodesic distances are com- 
puted using Floyd's algorithm for a directed graph JTsj . 
As might be expected from the absence of self-similarity 
in Fig. 1, we do not observe dynamical scaling in the 
above correlators over the time scales of our simulation. 
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Fig. 3 Growth laws extracted from the DTMC with 
Nf = 15A'^ (the straight line fits are aids to the eye), (a) 
Re ~ E'^^ crosses over to the usual t^^^ growth, while (b) 
Rhh exhibits a late time t^^^ behaviour. 

The length scales Re, Rh and Rhh (averaged over 16 
runs) are displayed in Figs. 3(a) and 3(b) for Nf = 15iV 
(i?0^ behaves the same as Re)- We identify two dynam- 
ical regimes associated with distinct growth processes. 
Regime I heralds the initiation of the strip instability to 
the formation of buds and the length scales, Re, Rh 
and Rhh grow as t^^* This is followed by regime 
II where the buds of typical size Rh move towards each 
other and subsequently coalesce. In the absence of any 
hydrodynamic drag due to the solvent. Re ~ t^^^. This 
dependence will be unchanged by inclusion of a viscous 
drag due to the solvent (though the prefactor would be 
altered). The buds carry both cj) and curvature H, and so 
equating the energy densities ~ i?^^ and Eh ^ R~hh^ 
leads to Rhh ^ Re^ ~ ^^^^ (Fig. 3(b)). 



It would be interesting to carry out a systematic ex- 
perimental study of dynamical morphology changes and 
shape instabilities in mixed artificial membranes. Of- 
course the choice of the lipid species is crucial; it is 
important to ensure that Tc > Tm, the main transi- 
tion temperature below which the lipid freezes into a gel 
state (Tm « 23°C for DMPC). This may be achieved 
by choosing lipids with short chains or with a lot of un- 
saturated bonds within the chain p9| . The dynamical 
buckling phenomenon reported above could be observed 
by inducing photochemical (UV) crosslinking of the mi- 
nority lipids after they have clustered. The present study 
might be of relevance to biological membranes as well, in 
particular to the dynamics of 'budding' of coated vesicles 
and the 'patching' of membrane proteins and glycolipids 
& 
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